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Abstract 

.^ . The usual power function error estimates do not capture the true 



order of uniform accuracy for thin plate spline interpolation to smooth 
data functions in one variable. In this paper we propose a new type of 
i-Q ■ power function and we show, through numerical experiments, that the 

C^ . error estimate based upon it does match the expected order. We also 

study the relationship between the new power function and the Peano 
kernel for univariate thin plate spline interpolation. 



1 Introduction 

For each 7 > 0, define the basis function 0^ : ]R — )■ ]R by 

xl^, if7 2IN, 



■^^^^^ 1 |x|^log|2;|, if7E2]N. 



Let m^ = [7/2J be the integer part of 7/2. For any integer n > ra^ and 
any set of values {fo, ■ ■ ■ ■, fn} of a target function / prescribed at the set of 
C^ i equi-spaced knots {0, /i, 2/i, . . . , 1}, where h = 1/n, Micchelli's theory [7| of 

conditionally positive definite radial basis functions guarantees the existence 
of a unique function Sh^^ of the form 

n "»7 

■5/1,7 (^) ~ z^ 0'k4>-y (^ ~ hk) + y^ bix , X G K, (1-1) 

fc=o 1=0 

that satisfies the interpolation conditions 

Sh,j{hi) = fi, i = 0, l,...,n, (1.2) 
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as well as the 'side' conditions 

n 

'^akk^ = 0, l = 0,l,...,m^. (1.3) 

fc=0 

In this paper we focus on the special case 7 = 2 corresponding to the thin plate 
spline (TPS) basis function 02 {x) = |xp log |x[ and we investigate the rate at 
which the interpolant (jl.ip converges to / uniformly over [0,1], as /i — )• 0. If 
/ has a Lipschitz continuous third derivative on [0, 1], Bejancu [3j proved that 
the uniform error over a fixed compact subset of (0, 1) inherits the maximal 
convergence rate 0{h^) obtained by Powell |9j and Buhmann [5] for 'cardinal' 
TPS interpolation on the infinite grid /iZ. Due to boundary effects, however, 
the uniform norm of the error over the full interval [0, 1] decays at the much 
slower rate 0(/i^'^), as illustrated numerically in P^HH [5]. 

The usual method for error estimation in radial basis function interpolation, 
reviewed in the next section, delivers bounds of the form 

\f{x) - Sh,-y (x) I < Cf^^Vh,-y{x), X E [0, 1], 

where Vh,'y is the so-called 'power function' associated with c/)^, while / belongs 
to the 'native space' generated by (p-y [ElIlT]. It is well known that theoretical 
convergence rates based upon bounding T'h,y{x) uniformly for x € [0, 1] do not 
match the actual rates of decay of the error achieved in numerical experiments 
if / has sufficiently many continuous derivatives. This discrepancy was first 
observed by Powell |10j for the bivariate TPS interpolant. 

For 7 = 2, in section 3, we obtain a new error bound which employs a 'mixed 
power function' M.h,fM defined by means of the basis functions (j)2 and 0^, for 
/U G (0, 4). We then perform a numerical study of max^-gjo,!] ■M.h,f_iix) as /i — )■ 0, 
which shows that, for ^u G [3, 4), the mixed power function decays like a constant 
multiple of /i^'^. This matches exactly the previously known numerical order 
of uniform convergence of the error / — Sh,2 on [0,1], for sufficiently smooth 
target functions /. In section 4 we prove that, for // = 3 and x G [0,1], the 
mixed power function value Aih,3ix) is, up to a constant factor, the L^-norm of 
the Peano kernel of the error functional at x. Moreover, we provide numerical 
evidence that the smaller L^-norm of this Peano kernel does not in fact decay 
faster than the mixed power function when measured uniformly over [0,1]. 
It is hoped that these results and the conjectures formulated in the paper 
will motivate future work to establish theoretically the uniform convergence 
order 0(/i'^") for univariate TPS interpolation to sufficiently smooth target 
functions. 



2 Error estimates via the standard power function 

In this section we review the power function technique to obtain error estimates 
for univariate interpolation with the radial basis function (p^ . A key role in this 
technique is played by the generalized or distributional Fourier transform of 

Lemma 2.1. [ISL section 8.3] For each 7 > 0, the generalized Fourier 
transform of (f>^ satisfies 

A, 

k|l+7' 



'P"f(i)=unh^ tGlR\{0}, 



for some constant A^ such that (— 1)"*^"^^^^ > 0. 

2.1 The standard power function 

As above, let m^ = [7/2J , n > m^, and h = 1/n. For each x G IR which is not 
in the knot-set {0, h, . . . ,1}, Micchelli's theory implies that the quadratic form 

(n n n \ 

Y^Y,^.Vf,<p^{hj-hk)-2Y,Vjct>.,{x-hj)\ , (2.1) 
j=0 fc=0 j=0 / 

is strictly positive whenever the non-zero vector v = [vq, . . . ,Vn) G IR"^^ 
satisfies 

n 

x^ = ^Vj{hjy, l = 0,l,...,m^. (2.2) 

j=0 



Further, for each j € {0,1,..., n}, let lyf^ be the unique function of the type 
(ll.ip - (ll.3p which satisfies the Lagrange interpolation conditions 



i^^liih) = 5ij, i = 0,l,...,n. 
Then we have the Lagrange representation formula for the interpolant (II. ip : 

n 

Sh,^{x) = Y,f/^l{x), xeIR, (2.3) 

j=o 

as well as the reproduction formula 

n 

^' = E(^-?')'4>(^)' ^^^' ^ = 0,l,...,m^. (2.4) 



Proposition 2.2. |17j With the above notations, for each x € IR, the vector 

T 



(^S(^)>--->^S(^)) ^^ 



n+1 



has the property that it minimizes the quadratic form (j2.ip among all non-zero 
vectors v G ]R"+i that satisfy (p:2]) . 

The minimum value of tiie quadratic form Q^,„ defines the square of the so- 
called 'power function' Vh,-y : IR — )• [0, oo), namely 

P^,^(x) := g^,„(v,). (2.5) 

Note that Vh,'y{hj) = 0, Vj € {0, 1, ... , n}. 
Proposition 2.3. |T7] For each x G IR, let 

n 

Q.,^{t):=e'^'-Y,^^-l{x)e'^'\ *eIR- (2.6) 

i=o 

T/ien we have the absolutely convergent integral representation 

^l,W = ^X^|#^*. (2.7) 

2.2 Error estimates 

For each 7 > 0, let k^ = [7/2 + 1] be the least integer that is greater than 
or equal to 7/2 + 1. In order to obtain error bounds for any target function 
/ G C"^t[0, 1], we construct an extension /* : IR ^- IR of / as follows (cf. [3j). 
By the Whitney extension theorem |16] . there exists / € C^t(IR) such that 
f{x) = f{x), for X G [0, 1]. Let u be an infinitely differentiable cut-off function 
which satisfies u^x) = 1 for x G [0,1] and z/(x) = for sufficiently large |x|, 
and set 

f*{x) := v{x)J{x), Vx G IR. 

Clearly /* G ^^(IR) is compactly supported and coincides with / on [0,1]. 
Furthermore, its Fourier transform /*, defined as the continuous function 

T*{t)--= I e-*^*r(x)dx. 



satisfies 



t^-f* (t) 



iP 



\(«7) 



(t) 



< 



if 



*\(k-,) 



Li(IR) 



(2.8) 



for any t 7^ 0. In particular, /* is integrable over IR, so /* can be recovered 
via the Fourier inversion formula 



fix) 



1 

2^ 



e*^7* (t) dt, xeM. 



(2.9) 



R 



Next, let fj := f {hj) in ([O]) for j = 0, . . . , n. Then ^^ and ^^ imply the 
error representation 



/ (x) - Sft,^ {x) = ^ f T* it) Qxa (t) dt, X G [0, 1] 



(2.10) 



where Qx,-) is given by (j2.6p . Moreover, as a consequence of (|2.8p and the 
definition of k^, we have 



C/,7 - 



|?(t)|2|t|l+^dt} 



1/2 



< 00, 



■2vr|yl^l j]R 

i.e., /* belongs to the so-called 'native space' generated by (f)^. Using (j2.7p and 
the Cauchy-Schwarz inequality in (j2.10p . we obtain the error bound 



1/ [x] - Sh,-y (x)| < Cf^^Vh,^ [x), X G [0, 1]. 



(2.111 



Further, Wu and Schaback [17] showed that the variational characterization of 
the power function given in Proposition 12.21 implies 



max Vh'yix) < B-yh^' , as /i ^ 0, 

xg[0,l] 



(2.12) 



for a constant B^ independent of h. On the other hand, Schaback and Wend- 
land [H] proved that the exponent 7/2 cannot be increased in the above bound. 
Thus, the power function technique leads to the maximal estimate 0{h^i'^) for 
the uniform norm of the error over [0, 1]. 



3 A mixed power function for univariate TPS 

In this section we focus on the TPS basis function (t)2{x) = |xplog|x|, i.e. 
7 = 2. According to ()2.12p . in this case the power function satisfies 



max Vh 2{x) = O (h) , as /i ^ 0. 

a;G[0,l] 



However, numerical experiments [21 [6l [TT] suggest that the uniform error 

max |/(x)-Sft,2(a;)| (3.1) 

a:G[0,lJ 

is of the magnitude of /i^' ^ for a sufficiently smooth target function /. 

To address this discrepancy, we start from the integral representation (|2.7p : 

^U-) = ^ / ^#^rf*- (3-2) 



Note that expression ()2.6p satisfies 



0(|tn, ast^O, 



|e.,2(t)|= ;:y W' "- -' (3.3) 

I C (1) , as \t\ — )■ oo. 



for each fixed x and /i. Indeed, since m2 = 1, by ()2.4p we have 

1 = ^ £J5 (^) a^d ^ = IZ ^•^' ■ 45 (^) ' for ^ ^ ^- (3-4) 

i=o ' i=o 

This fact, together with the series expansion of the exponential, provides the 
bound (j3.3p for t — )• 0. The bound for |t| — )• cx) follows from the triangle 
inequality. 

As a consequence of (j3.3p . the integral (j3.2p is still well defined if |tp is replaced 
by 1*1"'^"'''^, for any /x G (0,4), [x ^ 2. We may thus define the mixed power 
function Mh,tJ. : IR — > [0, oo) whose square is given by 

Ml^ (x) := ^ ^ %^^*' ^ ^ ^' ^ ^ (0' 4) • (3-5) 

Under this integral, the Lagrange functions entering in expression (j2.6p of 
0x,2 (and generated by the TPS basis function (j)2) are combined with the 
generalized Fourier transform of the basis function c/)^ (cf. Lemma 12. ip . 

We now let n = \^/2 + 1] , / S C"^ [0, 1] and, as in subsection 2.2, consider the 
compactly supported extension /* € C^ (IR) of / to the whole real axis. Then 
the error analysis of subsection 2.2 recast in terms of the mixed power function 
implies 

1/ {x) - Sh,2 {x)\ < Cf^^Mh,^ {x), xe [0, 1], fie (0, 4) , (3.6) 



where 

f 1 /■ -— - 2 ^^ -,1/2 

This shows that, for a fixed ^ G (0,4), esthnates of the decay of the mixed power 
function Mh,fi as /i — )• wih dehver error estimates for TPS interpolation. 
Therefore we state the fohowing problem: 

Problem 3.1. Given /i G (0,4), fJ- ^ 2, does there exist an algebraic decay rate 
of the mixed power function uniformly over [0, 1], i.e., a largest value a^ > 
such that 

max A^/,,^(x) = 0(/i"''), as /i ^ 0? (3.8) 

xG[0,1] 

Before embarking on a numerical answer to this problem, a few remarks are in 
order. Firstly, note that, due to (13. 7p . the above target function / e C^ [0, 1] 
has its compactly supported extension /* in the native space generated by the 
basis function (pfj_, rather than the native space generated by the TPS basis 
function (j)2 as in the standard estimate (j2.1ip for 7 = 2. Thus, for a given 
// G (0,4) {fi y^ 2), the resulting mixed power function error bound (|3.6p is 
precisely what we would expect if we measured the TPS interpolation error for 
target functions in the native space of (/>^; this approach is investigated in [12] 
in the context of approximation rather than interpolation. In particular, the 
mixed power function bound (j3.6p applies to the smooth (C°°) target functions 
employed in the numerical experiments of [21 El [H] . 

Secondly, recall the two equivalent expressions for the standard power function: 
the direct form (12. 5p and the integral representation ()2.7p . Letting m = [/x/2j , 
an application of Theorem 3 from [T7] shows that the mixed power function 
can also be expressed as 

(n n 
j=0 k=0 

(3.9) 

j=o J 

Thirdly, note that ()3.4p implies that the TPS Lagrange functions £■ ^ satisfy 
constraint (12. 2p of the variational problem from Proposition 12. 21 with /i in place 
of 7. However, the solution to that problem is provided by the values of the 





/x = l/3 


/x = 2/3 


/" = ! 


h-' 


j^{max) 


ah,fi 


M^r' 


ah,fi 


M^r^ 


ah,fi 


128 


4.774E-01 


0.167 


1.768E-01 


0.333 


6.342E-02 


0.500 


256 


4.253E-01 


0.167 


1.404E-01 


0.333 


4.485E-02 


0.500 


512 


3.789E-01 


0.167 


1.114E-01 


0.333 


3.171E-02 


0.500 


1024 


3.376E-01 


0.167 


8.842E-02 


0.333 


2.242E-02 


0.500 


2048 


3.007E-01 


0.167 


7.018E-02 


0.333 


1.586E-02 


0.500 


Cf, 


1.072 


0.8912 


0.7175 



Table 1: Decay of the mixed power function for fi € (0, 1] 



Lagrange functions generated by interpolation with the basis function 0^. As 
a result, the bounding technique [ITj that leads to the estimate (j2.12p for the 
standard power function cannot be applied to obtain estimates on Mh,^- 

We now turn to a numerical investigation of the behaviour of the mixed power 
function. For a fixed parameter ^u G (0,4), /x 7^ 2, we compute an approxima- 
tion -M-f^i^^ of the left-hand side of (j3.8p for h = 1/n, starting from n = 128 
and proceeding as follows: 



1. For the current mesh-size h and each jG{0,l,...,n}, express the TPS 
Lagrange function £■ ^ in the form (jl.ip and compute its coefficients by 
solving the system (|l.2p - (jl.3p . where fi = 6ij, i G {0, 1, . . . , n}. 

2. Use (j3.9p to evaluate the mixed power function at the set of mid-points 
^evai,h = {h/2, 3/1/2, . . . , 1 — h/2} and determine its maximum value 



M 



(max) 
h,ij, 



max{Mh,^iix) : x G Xeval,h} 



3. Double n and repeat steps 1-2 as long as n < 2048. 



The results displayed in Tables [IHH show that, for each chosen /i, the values of 



Mi',"'^' satisfy 



M 



(max) 
h,ii 



c^/l"'•■^ 



where c^ and at,^ are also included in the tables. On the basis of these nu- 
merical results, we are led to the following conjecture. 





fi = 4/3 


;U = 5/3 


h-' 




Oih,fi 


j^(max) 


ah,,! 


128 

256 

512 

1024 

2048 


2.143E-02 
1.350E-02 
8.503E-03 
5.356E-03 
3.374E-03 


0.667 
0.667 
0.667 
0.667 
0.667 


6.258E-03 
3.512E-03 
1.971E-03 
1.106E-03 

6.208E-04 


0.833 
0.833 
0.833 
0.833 
0.833 


Cf, 


0.5442 


0.3568 



Table 2: Decay of the mixed power function for /x G (1,2) 





/i = 7/3 


;U = 8/3 


H = 3 


h-' 


Mtr 


a/i,^ 


M^r^ 


a/i,^ 


MiT 


ah,tM 


128 


1.061E-03 


1.167 


6.221E-04 


1.334 


3.327E-04 


1.507 


256 


4.727E-04 


1.167 


2.473E-04 


1.334 


1.196E-04 


1.503 


512 


2.106E-04 


1.167 


9.828E-05 


1.333 


4.296E-05 


1.500 


1024 


9.381E-05 


1.167 


3.905E-05 


1.333 


1.543E-05 


1.498 


2048 


4.179E-05 


1.167 


1.551E-05 


1.333 


5.534E-06 


1.496 


c^ 


0.3049 


0.4024 


0.4975 



Table 3: Decay of the mixed power function for ^ € (2, 3] 





/i = 10/3 


/i = 11/3 


h-^ 




"/i,At 


j^.(max} 


Oih,fj. 


128 


2.032E-04 


1.491 


1.661E-04 


1.497 


256 


6.995E-05 


1.497 


5.808E-05 


1.499 


512 


2.419E-05 


1.501 


2.039E-05 


1.500 


1024 


8.402E-06 


1.503 


7.182E-06 


1.501 


2048 


2.919E-06 


1.505 


2.523E-06 


1.502 


Cm 


0.2814 


0.2368 



Table 4: Decay of the mixed power function for /i G (3, 4) 



Conjecture 3.2. The mixed power function satisfies the estimate ()3.8p with 
the algebraic decay rate 



f, /or /xG (0,3) \ {2}, 
I, /or/xG[3,4). 



4 The mixed power function for fi = 3 

Note that if Coniecture 13.21 can be established for a particular value ;U G [3,4), 
then the mixed power function bound ()3.6p implies a new and improved error 
estimate for thin plate spline interpolation on the unit interval, namely, for any 

/gC3[o,i], 

max \f(x)-Sh2ix)\ = 0(h^/A , as /i ^ 0. (4.1) 

xG[0,l] ' V / 

This would provide a theoretical explanation of the numerical results reported 

in [21 El [n]. 

In this section, we investigate Conjecture 13.21 for the special case fj, = 3. By 
(13. 5p and ()3.9p . the square of the mixed power function A4h,3, combining TPS 
Lagrange functions with the cubic basis function ^3, is given by 

MUr, . ^ f ^^.< 

n n n V^-^j 

Figure [1] illustrates the decay of Mh,2, as /i ^- 0. It can be confirmed nu- 
merically that the decay rate 0{h^''^) of M.h,3 suggested by Table [3] applies 
uniformly on [0, 1], i.e. all peaks of the plot decay at this rate. 

We now relate the mixed power function Aih,3 to a classical error analysis 
method, namely the Peano kernel representation. Let 

n 
j=0 

For each x G [0, 1], Eh^x is a continuous linear functional on C[0, 1] with the 
usual max norm, and (|3.4p implies that the linear polynomials are in the null 
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1 /h = 16 
1 / h = 32 
1 / h = 64 



\MMMmmMmm 



Figure 1: Plot of Mh,3 (x) for h'^ = 16, 32 and 64. 



space of Efi^x- Then, for any / with an absolutely continuous first derivative 
on [0,1], Peano's theorem [U p. 271] implies 

f{x) - Sh,2{x) = [ Kh,Au)f"{u) du, Vx G [0, 1], (4.3) 

Jo 

where K^^^ is the 'Peano kernel' given by 

n 

KhAu) ■■= {x - u)+ - Y,i^ll{x){hj - u)+, n G R. 

j=o 

Proposition 4.1. For each x G [0, 1], the mixed power function value 
■M.h,3 {x) is a constant multiple of the L'^[0, l]-norm of the Peano kernel Kh^x- 

Proof. The reproduction property ()3.4p implies that -K'h.x is compactly sup- 
ported on [0, 1], and that 



Kh,x{u) = T,\\x-u\- 'S2^jhi^)\^J -u\\ , 



n G IR. 



Then, by Lemma 12. H the Fourier transform of the kernel K^^^ is the analytic 
and square integrable function 

11 



where Qx.2 is defined by (j2.6p . Therefore, using the first hne of (|4.2p . the 
Parseval-Plancherel formula, and the compact support of Kj^x-, we deduce 

^1 Jtr 

which is the required conclusion. D 

As a consequence, we obtain an alternative way of bounding the error f — Sh,2 in 
terms of the mixed power function A4h,3, by using Cauchy-Schwarz directly in 
the right-hand side of the Peano formula (14. 3[) . The resulting bound applies to 
any / with an absolutely continuous first derivative on [0, 1] and /" G i^[0, 1]. 
This represents an improvement over (|3.6p - (l3.7p . which required / € C'^[0, 1] 
for /x = 3. 

Finally, a related question of interest is whether a sharper uniform error bound 
can be obtained from (j4.3p via Holder's inequality 

\f{x) - Sh,2{x)\ < Bh{x) ||/"|L..[o,i] ' ^ ^ [0' 1]' 
where/" E L~[0, 1] and 

Bh{x) := / \Kx^hiu)\du. 
Jo 

We note that this technique was used by Atkinson [T] in the late 1960s to inves- 
tigate the error behavior of natural cubic spline interpolant near the endpoints 
of the unit interval; see also Schaback |13] for a treatment that is closer to our 
presentation. In the case of the TPS interpolant, a numerical answer to the 
question is provided in Table O whose entries satisfy 

Bh(-j =0.05059 h'^'', and Bh (^—] = O.U955h'"^ , 

i.e., Bh decays approximately with the rate 0{h^''^) near the endpoints of [0, 1] 
and this rate improves to 0{h'^) near the midpoint. Also, Figure [2] shows that 
the extreme peak value is well approximated by Bh (2)1 while all of the lower 

12 



peaks decay at the faster rate. Therefore estimating the L^-norni of the Peano 
kernel leads to the same rate of decay 0(/i^'^) of the uniform error (j3.ip as 
that predicted in (|4.ip by the mixed power function A4h,3- 

We conclude the paper by remarking that any theoretical proof of the uniform 
decay rate 0{h^''^) of Aih,3{x) or Bh{x) for x G [0, 1] will have to rely on specific 

(2) 

properties of the TPS Lagrange functions £■ ^, j G {0, 1, . . . , n}. A potentially 
useful such property is the special case of [H Theorem 3.1] stating that the 
Lebesgue-type constant 



max 

xe[o. 



n 



j=0 



admits an upper bound independently of the mesh-size h. It remains an open 
question whether this or other properties of the TPS Lagrange functions can 
lead to further progress on the above conjectures. 



'- 1 /ii = 16 

- - 1 / h = 32 

1 / h = 64 



Sm^MNSN^I'fhl'MNW^^ 



Figure 2: Plot of Bh{x) for h'^ = 16, 32 and 64 
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/l-l 


^41) 


h 


Bh{^) 


o-fe 


64 


1.024E-04 


1.491 


3.633E-05 


2.001 


128 


3.533E-05 


1.498 


9.098E-06 


2.001 


256 


1.228E-05 


1.501 


2.293E-06 


1.999 


512 


4.289E-06 


1.503 


5.694E-07 


2.000 


1024 


1.502E-06 


1.504 


1.434E-07 


1.999 



Table 5: Decay of the L^-norm of the Peano kernel. 
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